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Fusion research on magnetic confinement is 
confronted with a severe problem concerning the 
electron densities to be used in fusion de- 
vices. Indeed, high densities are mandatory for 
obtaining large efficiencies, whereas it is empir- 
ically found that catastrophic disruptive events 
occur for densities exceeding a maximal one n*^. 
On the other hand, despite the large theoretical 
work "there is no widely accepted, first principles 
model for the density limit" (see [1], abstract). 
Here, we propose a simple microscopic model of 
a magnetized plasma suited for a tokamak, for 
which the existence of a density limit is proven. 
This property turns out to be a general collective 
feature of electrodynamics of point charges, which 
is lost in the continuum approximation. The law 
we find is 
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where /io is the vacuum permeability, c the speed 
of light, 77ie the electron mass, and B the magnetic 
field. As shown in Fig. 1, the theoretical limit 
(big circles) is in rather good agreement with the 
empirical data, actually a surprisingly good one 
for a model based on first principles, with no ad- 
justable parameter. 

The way in which law ([l} was established is an inter- 
esting example of an encounter between fundamental and 
applied research. Indeed three of the present authors are 
involved since some time in studies of a general character 
concerning the microscopic electrodynamics of systems of 
point particles (see [3] and [s^l), in which both the mu- 
tual retarded forces, and the radiation reaction force of 
Abraham Lorentz and Dirac [4| (see also [5| and [(|) are 
taken into account. One of the results obtained is the 
proof of an identity conceived by Wheeler and Feynman 
Q, and an appreciation of the role the latter plays in 
allowing for the very existence of a dispersion relation. 
In particular, some examples of dispersion relations were 
given (see [3], Fig. 1), which exhibit, as the matter den- 
sity is increased, a bifurcation of a topological character, 
entailing an instability. But the physical relevance of this 
fact was not emphasized. Such a density controlled bifur- 
cation impressed instead very much those of the present 
authors who deal with plasma physics, who suggested it 
may be relevant for fusion plasmas. To this end, the sim- 
plest possible model was formulated, that should capture, 
within the frame of the foundational works mentioned, 
the essential physics of a magnetized plasma, confined in 
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FIG. 1: Greenwald plot for three tokamak devices (data ex- 
tracted from [Jl ) . Dotted line is the empirical Greenwald limit 
@ . Big circles are the theoretical predictions ^ for the same 
devices. 



a tokamak configuration (the most studied one for fusion 
plasmas). The model is presented here, together with the 
deduction of law ([1]) . Preliminarily, the main evidence for 
the existence of a density limit in tokamaks is recalled, 
and it is discussed how well does law ((1} fit the data. 

For the purposes of the present paper, all is needed to 
know about tokamaks is essentially that they are toroidal 
devices in which the confining magnetic field B is the 
vector sum of a strong toroidal field Bt produced by a 
set of coils wound around a torus, and of a much smaller 
poloidal field Bp generated by a toroidal plasma current 
Ip. A few more details will be mentioned later. 

The empirical data which show the existence of a den- 
sity threshold beyond which tokamaks cannot operate, 
were collected by Greenwald fll] in a classical figure, the 
data of which are reported here in Fig.l. In the figure, the 
electron densities ne at which three diff'erent tokamak de- 
vices could actually be operated are reported versus the 
so-called Greenwald parameter uq = Ipj'Kc? (where a is 
the minor radius of the torus), which is presumed to be 
the relevant control parameter. Indeed, the law proposed 
by Greenwald for the maximal density (dotted line in the 
figure) is 



rig = ac Ip/i^a 



(2) 



where ac is a constant with suitable dimensions, such 
that ac — 1 in the units indicated in the figure. The 
theoretical predictions given by ([TJ are also reported as 
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big circles. 

Thus there naturally arises the question, how is it pos- 
sible that two analytically different predictions, ([Ij and 
(HI), happen to agree with each other, at least in a few 
definite cases. The reason is that one has 
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where the dimensional coefficient h (a length) is not a 
universal one, but depends parametrically on geometric 
factors and operative conditions characterising each ex- 
periment. These are the major and the minor radii R 
and a, the plasma elongation k (equal to 1 for circular 
plasmas, see Q pg. 277), and the edge safety factor q{a) 
defined below. 

This is seen as follows. One has Bp — jiQlp/ {2naK,^/^). 
Furthermore, in the approximation B Bt, one has 
B « {Bt/ Bp)Bp, while the operative parameter Bt/Bp is 
determined by the edge safety factor q{a) which, for the 
simplified case of an elliptical plasma (see [8] ) , is defined 
as q{a) = {a/R){Bt/Bp)(l + k^) /2. This gives relation 
The formula for h in the general case is also easily 
established. 

Now, magnetohydrodynamic stability requires q{a) > 
2, but the actual value at which each experimental data 
point of Fig. [T] was taken is not given in the literature. 
So we assumed q{a) = 4, which is a typical operational 
value, and this introduces an uncertainty in the theoret- 
ical points reported in Fig. [TJ 

Thus, we decided to look directly at the experimental 
values available in the literature, from which a definite 
estimate of B could be obtained, and this we did for 
conventional tokamaks more recent than those in Fig. 
HJ Such values are plotted in Fig. where they are 
compared to laws ((T} and (0). The agreement with law 
dl]) is perhaps a little better. 

This fact might have relevant implications for future 
tokamaks, as it implies a favorable density scaling for 
machines with large values of the product RB. For ex- 
ample, let us consider the international thermonuclear 
experimental reactor (ITER) '21'], which should operate 
at a toroidal field of 5.3 T. According to law ([1]) it would 
be able to operate at densities up to = 4.7x 10^°m~'^. 
This is a value more than three times larger than that 
expected according to Greenwald law ^ for a plasma 
current of 15 MA, which is the corresponding value of Ip 
for a q{a) = 3 scenario. 

It is worth mentioning that a B^ dependence of the 
density limit in the ALCATOR C experiment was noted 



in the past by Granetz [22|, although such a clear de- 
pendence was not observed on other experiments (see for 
example [13] for the DIII-D tokamak). Also, one should 
point out that tokamaks with very low aspect ratio R/a 
(spherical tokamaks) and reversed field pinches (RFP) 
|23| . for which the validity of the Greenwald scaling has 
been proposed , seem not to fit well into the proposed 
B^ scaling. This, perhaps, suggests that at low magnetic 
fields other effects, not considered in the present simple 
model, might come into play. 

We show now how law ([1} was obtained, in the frame 
of microscopic electrodynamics of point particles (see 
and Q), rather than of magnetohydrodynamics (see for 
example [2^, chapter 17), or of the mean field theories of 
the Vlasov approach. We describe the plasma as consti- 
tuted of point particles obeying Newton equations, with 
both the retarded electromagnetic interactions among all 
particles and the radiation reaction force taken into ac- 
count. We then concentrate on the role played by the 
gyration of the electrons around the magnetic field lines, 
and so ignore their motions along the field lines, and also 
ignore the electrostatic part of the problem. Finally we 
also limit ourselves to the extremely simplified case of a 
one-dimensional array. 

So we introduce the following model. Given a con- 
stant magnetic field B, which we take oriented along 
the z axis of a cartesian coordinate system, we con- 
strain each electron, say the n-th one, to move on a 
plane parallel to the {x^y) plane, so that its Zn coor- 
dinate is fixed. The simplest choice is to take z„ = n/, 
with n G Z, for a given positive step /. Each electron, 
say the n-th one, is subjected to the external magnetic 
field B, and also to the electromagnetic field created 
by the electrons themselves. Namely, the sum of the 
Lienard-Wiechert fields, which are determined as the re- 
tarded solutions of Maxwell equations having as sources 
the charge and the current densities of each other elec- 
tron m 7^ n, and the Abraham-Lorentz-Dirac radiation 
reaction force, due to the motion of the n-th electron 
itself. The latter is given, in the nonrelativistic approx- 
imation, by (2/3)e^/(47reoc^)xn, where e is the electron 
charge and x„ the position vector of the electron. 

We then perform the dipole approximation. Thus we 
neglect the magnetic field due to the m-th electron, and 
for the electric field created by it we take the well known 
expression for a dipole. Finally, we approximate the dis- 
tance between electrons n and m by rn,m = l\n- — 
The system of equations of motion defining the model is 
then 
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FIG. 2: Density limit values for various conventional tokamak machines us B in a and vs Ip/na^ in b. Dotted line in a is the 
proposed theoretical law ([T}, and in b is the empirical law For the various devices the points given are representative of 
disruptive events actually due to an increase of density or of an operative condition declared to be close to a density limit (see 
references P-[20|'). 



for n E Z, where eg is the vacuum permittivity, and 
LOc = eB/rUe the Larmor or cyclotron frequency of the 
electrons in the external magnetic field B. This is an 
infinite system of linear equations with delay, which is 
just a simple variant of the system considered in 

Our aim is now to investigate the stability properties 
of the system, as the control parameters rie — l/l^ and B 
(or equivalent ly ujc) are varied. Following a completely 
standard procedure (see for example [11]), we compute 
the normal modes of the system and determine the val- 
ues of the parameters for which the frequencies become 
complex. So we look for normal mode solutions with 
wavenumber k and angular frequency w, i.e., of the form 
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This leads to a linear system in the unknowns A^, Ay, 
from which the dispersion relation between uj and k is 
found by equating the determinant to zero. This gives 
two real equations in the two unknowns lo and fc, namely. 
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± — +pF(fc/,/a;/c) = 0, p^^. (6) 
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uj^ - G{kl,luj/c) = 
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Here, Up is the familiar plasma frequency defined by 

Wp = e'^/eomj^ = Ue e^/eoWe , (8) 

while F and G, as functions of the variables a = kl, 

j3 — luj / c, are defined by 



G(a,/3)=/33-<7(a,/3) , 



/32log(2|cos/3-cosa|) -/(a,/3) 



the functions / and g being the ones already introduced 
in [2:], namely. 



,cos{na — |n|/3) sin(na — \n 



-) 



.sin(na — \n\f3) ^ cos(na — |n|/3) . 



Some details concerning the summation of the series lead- 
ing to the term log (2 | cos /3— cos a |) entering the func- 
tion F are here omitted. 

Now, one meets here with a deep question of principle. 
Indeed, for fixed values of the parameters I and uJc one has 
two equations in two unknowns (w and k), and this would 
not allow for the existence of a dispersion relation, i.e., 
of a function uj — a;(fc) for a continuous range of values 
of k. However, the existence of a dispersion relation is 
guaranteed by the fact that equation ([7]) actually is an 
identity. In fact, this is a particular case of a general 
identity, conceived by Wheeler and Feynman Q and first 
proven in (see section 6) for a one-dimensional case 
and in Q for a three-dimensional one. 

So, the problem of obtaining the dispersion relation 
is reduced to solving ^ in the unknown uj — uj{k), 
in which p ~ uip/uj'^ plays the role of a parameter. In 
Fig. [3] the dispersion relations are shown for a cyclotron 
frequency Wc = 3.8 ■ 10^^ Hz, and for several values of 
the parameter p — ujp/uj'^ (or of the corresponding elec- 
tron density ne/nc, normalized to the Greenwald density 
no = Ip/ira^). 

The most important qualitative result is that normal 
modes are found to exist (for all k) only below a critical 
value of p, i.e., below a certain threshold of plasma den- 
sity. Indeed, starting up from low densities, at a certain 
critical density a bifurcation is seen to occur, character- 
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FIG. 3: The dispersion curves, solutions of equation © (with 
the plus sign), in the plane {kl, uj/uc) for oJc = 3.8 x lO'^^ Hz, 
and for several values of the parameter p = LUp /uj^ (or equiva- 
lently of the electron density rie, normalized to the empirical 
no limit). 

ized by the fact that the curves no more intersect the 
vertical axis kl = n. This means that for values of k just 
below tt/I equation ^ does not admit a real solution, so 
that the corresponding frequencies acquire an imaginary 
part, and the whole system becomes unstable. Numerical 
computations not reported here show that the character- 
istic time of the instability is of the order of 1-k jiOc and 
that above the critical density the Wheeler and Feynman 
identity is no more satisfied. 

Notice that this phenomenon of the existence of a max- 
imal allowed density is obviously lost if one introduces the 
continuum approximation, i.e., is a characteristic feature 
of the discrete structure of matter. Indeed, following Q, 
the continuum approximation corresponds to deal with 
wavelengths much larger than the step i.e., to assume 
k tt/^, whereas the existence of a density limit depends 
on the behavior of the system for kl ~ tt. 

We have now to determine the bifurcation value of the 
parameter p. As the bifurcation occurs for kl = n and 



for values of w/wc < 1, i.e., for luj/c ~ 0, one can just 
limit oneself to study equation ^ for a fixed value of the 
function F, namely F{Tr, 0), so that one is simply reduced 
to deal with an algebraic equation of second degree. One 
computes F(7r, 0) ~ 0.14, and so real values of uj are 
found to exist only for p < 1.74. This, together with the 
definition of p in ([6]) and eo^o = gives law 

Notice that law H]) has the same form of the Brillouin 
limit (27l |. which is known to apply to the case of non- 
neutral plasmas [l^ . The main difference with respect to 
our procedure is that in the case of the Brillouin limit the 
electric field acting on each electron is introduced within 
a mean field approach, whereas here it is computed in 
the frame of a many-body microscopic theory. Corre- 
spondingly, we find that the instability involves normal 
modes with wavelengths of the order of the mean elec- 
tron distance, so that it escapes a mean field approach. 
In particular, such an instability is found to occur in neu- 
tral plasmas, for which the mean charge density vanishes, 
and the Brillouin approach cannot be used. 

A final comment concerns the possibility of dealing 
with the other main magnetic configuration studied for 
the confinement of fusion-relevant plasmas, i.e., the Stel- 
larator 0|. The present model does not directly apply. 
Indeed, in the Stellarator a large amount of power is typ- 
ically transferred to the electrons through electron cy- 
clotron resonance heating (ECRH), and this requires to 
add in our model a forcing term. 

In conclusion, through an extremely simplified model 
of a magnetized plasma suited for a tokamak, based on 
first principles, we have proved the existence of a density 
limit, beyond which the system becomes unstable. The 
law thus found differs from the usually accepted one, and 
this fact might have relevant implications for future toka- 
maks. 
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